Interaction QTL x Year
Interaction QTL x Year
This file aims to provide more details concerning the QTL x Year interaction of the resistance of Durum Wheat to WSSMV. We are going to study ELISA
library(QTLRel)## R/QTLRel is loaded
library(xtable)Introduction
Can we do exactly the same thing with lm instead of QTL-Rel? The only difference is that we do not take the population structure into account.
Load data
I load the genotyping matrix first.
genotype<-read.table("/Users/holtz/Dropbox/Publi_Mosaique/DATA/DATA/GROUPED/genotypage.csv", sep = ";" , header = F, na.strings = "-")
genotype=as.matrix(genotype)
colnames(genotype)=genotype[1,]
genotype=as.data.frame(genotype[-1 , ])
names(genotype)[1]<-"geno"
print("--- Your genotyping matrix looks correct. Dimension of the matrix are :")## [1] "--- Your genotyping matrix looks correct. Dimension of the matrix are :"
print(dim(genotype))## [1] 348 7342
# I copy this matrix 2 times, since I read 2012 and 2015 together.
rownames(genotype)=genotype[,1]
genotype=genotype[,-1]
a=genotype ; rownames(a)=paste(rownames(a),"2012",sep="_")
b=genotype ; rownames(b)=paste(rownames(b),"2015",sep="_")
genotype=rbind(a,b)Then the genetic map:
map <- read.table("/Users/holtz/Dropbox/Publi_Mosaique/DATA/DATA/genetic_map.txt" , header=T , dec = ".", na.strings = "-" , check.names=F)
colnames(map) <- c("LG", "marqueur", "Distance","group_physique","Posi_physique")
rownames(map) <- map$marqueur
map$LG <- as.factor(map$LG)
print("--- Your genetic map looks correct. Dimension of the map are :")## [1] "--- Your genetic map looks correct. Dimension of the map are :"
print(dim(map))## [1] 8568 5
map=map[ , c(2,1,3,5)]
colnames(map)=c("snp","chr", "dist", "phyPos")And finally the phenotyping matrix
BLUP<-read.table("/Users/holtz/Dropbox/Publi_Mosaique/DATA/DATA/GROUPED/phenotypage.csv", header = TRUE, sep=";")
colnames(BLUP)[1]="geno"
print("--- Your Phenotyping matrix looks correct. Dimension of the matrix are :")## [1] "--- Your Phenotyping matrix looks correct. Dimension of the matrix are :"
print(dim(BLUP))## [1] 396 6
# Fichier de phénotypage modifié, il va falloir mettre la Elisa de 2012 et 2015 ensemble, avec une colonne année.
a=BLUP[, c(1,3)] ; a$year="2012" ; colnames(a)=c("geno", "Elisa_blup_AR1","year")
#a[,2]=a[,2]/sqrt( mean(c(0.66,0.83)) )
#a[,2]=ifelse(substr(a$geno, 1,2)=="TT", a[,2]/sqrt(0.66) , a[,2]/sqrt(0.83))
b=BLUP[, c(1,5)] ; b$year="2015" ; colnames(b)=c("geno", "Elisa_blup_AR1", "year")
#b[,2]=b[,2]/sqrt( mean(c(1.2,1.21)) )
#b[,2]=ifelse(substr(b$geno, 1,2)=="TT", b[,2]/sqrt(1.2) , a[,2]/sqrt(1.21))
BLUP=rbind(a,b)
rownames(BLUP)=paste( BLUP[,1], BLUP$year,sep="_")
BLUP=BLUP[,-1]
# Note: On peut garder les blups tels quels / ou les pondéré par la variance génet de chaque année moyennée sur les 2 pops / ou par la variance génet de chaque année et chaque pop.Prepare data
We need to have genotype and phenotype in the same order.
And I add a “pop” column in the phenotyping matrix:
Y=na.omit(BLUP)
Y=Y[which(rownames(Y)%in%rownames(genotype)) , ]
Y$pop=substr(rownames(Y),1,2)
genotype=genotype[which(rownames(genotype)%in%rownames(Y)) , ]
genotype=genotype[ match(rownames(Y),rownames(genotype)) , ]test1: DS - 2012
QTL detection without QTL rel for DS only in 2012 only
# Initialize the result table
result_lm_2012_DS=data.frame(matrix(0,0,3))
colnames(result_lm_2012_DS)=c("marker","pval-marker","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$year=="2012") , ]
don=don[which(don$pop=="TT") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_DS[num,1:3]=c(i, res$`Pr(>F)`[1],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_DS[,2:3]=apply(result_lm_2012_DS[,2:3] , 2 , as.numeric)
dim(result_lm_2012_DS)## [1] 3544 3
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODs with the genetic map
result_lm_2012_DS=merge(map,result_lm_2012_DS, by.x=1 , by.y=1, all.y=T)
result_lm_2012_DS=result_lm_2012_DS[order(result_lm_2012_DS$chr, result_lm_2012_DS$dist) , ]
# And plot it
plot(-log10(result_lm_2012_DS$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_DS$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_DS))
num=aggregate(num, by=list(result_lm_2012_DS$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test2: DS - 2015
QTL detection without QTL rel for DS only in 2012 only
# Initialize the result table
result_lm_2015_DS=data.frame(matrix(0,0,3))
colnames(result_lm_2015_DS)=c("marker","pval-marker","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$year=="2015") , ]
don=don[which(don$pop=="TT") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele)
res=anova(model)
# Add result to the 'result' file
result_lm_2015_DS[num,1:3]=c(i, res$`Pr(>F)`[1],summary(model)$r.squared)
}
}
# wrong class
result_lm_2015_DS[,2:3]=apply(result_lm_2015_DS[,2:3] , 2 , as.numeric)
dim(result_lm_2015_DS)## [1] 3544 3
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODs with the genetic map
result_lm_2015_DS=merge(map,result_lm_2015_DS, by.x=1 , by.y=1, all.y=T)
result_lm_2015_DS=result_lm_2015_DS[order(result_lm_2015_DS$chr, result_lm_2015_DS$dist) , ]
# And plot it
plot(-log10(result_lm_2015_DS$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2015_DS$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2015_DS))
num=aggregate(num, by=list(result_lm_2015_DS$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test3: DL - 2012
QTL detection without QTL rel for DL only in 2012 only
# Initialize the result table
result_lm_2012_DL=data.frame(matrix(0,0,3))
colnames(result_lm_2012_DL)=c("marker","pval-marker","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$year=="2012") , ]
don=don[which(don$pop=="BX") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_DL[num,1:3]=c(i, res$`Pr(>F)`[1],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_DL[,2:3]=apply(result_lm_2012_DL[,2:3] , 2 , as.numeric)
result_lm_2012_DL=na.omit(result_lm_2012_DL)
dim(result_lm_2012_DL)## [1] 5851 3
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LOD with the genetic map
result_lm_2012_DL=merge(map,result_lm_2012_DL, by.x=1 , by.y=1, all.y=T)
result_lm_2012_DL=result_lm_2012_DL[order(result_lm_2012_DL$chr, result_lm_2012_DL$dist) , ]
# And plot it
plot(-log10(result_lm_2012_DL$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_DL$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_DL))
num=aggregate(num, by=list(result_lm_2012_DL$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test4: DL - 2015
QTL detection without QTL rel for DL only in 2015 only
# Initialize the result table
result_lm_2015_DL=data.frame(matrix(0,0,3))
colnames(result_lm_2015_DL)=c("marker","pval-marker","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$year=="2015") , ]
don=don[which(don$pop=="BX") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele)
res=anova(model)
# Add result to the 'result' file
result_lm_2015_DL[num,1:3]=c(i, res$`Pr(>F)`[1],summary(model)$r.squared)
}
}
# wrong class
result_lm_2015_DL[,2:3]=apply(result_lm_2015_DL[,2:3] , 2 , as.numeric)
result_lm_2015_DL=na.omit(result_lm_2015_DL)
dim(result_lm_2015_DL)## [1] 5851 3
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LOD with the genetic map
result_lm_2015_DL=merge(map,result_lm_2015_DL, by.x=1 , by.y=1, all.y=T)
result_lm_2015_DL=result_lm_2015_DL[order(result_lm_2015_DL$chr, result_lm_2015_DL$dist) , ]
# And plot it
plot(-log10(result_lm_2015_DL$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2015_DL$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2015_DL))
num=aggregate(num, by=list(result_lm_2015_DL$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])test5: DS AND DL - 2012
QTL detection without QTL rel for DS and DL in 2012 only
# Initialize the result table
result_lm_2012_DSDL=data.frame(matrix(0,0,4))
colnames(result_lm_2012_DSDL)=c("marker","pval-marker","pval-pop","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$year=="2012") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele + don$pop)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_DSDL[num,1:4]=c(i, res$`Pr(>F)`[1:2],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_DSDL[,2:4]=apply(result_lm_2012_DSDL[,2:4] , 2 , as.numeric)
dim(result_lm_2012_DSDL)## [1] 7341 4
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODs with the genetic map
result_lm_2012_DSDL=merge(map,result_lm_2012_DSDL, by.x=1 , by.y=1, all.y=T)
result_lm_2012_DSDL=result_lm_2012_DSDL[order(result_lm_2012_DSDL$chr, result_lm_2012_DSDL$dist) , ]
# And plot it
plot(-log10(result_lm_2012_DSDL$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_DSDL$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_DSDL))
num=aggregate(num, by=list(result_lm_2012_DSDL$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test6: DS AND DL - 2015
QTL detection without QTL rel for DS and DL in 2015 only
# Initialize the result table
result_lm_2015_DSDL=data.frame(matrix(0,0,4))
colnames(result_lm_2015_DSDL)=c("marker","pval-marker","pval-pop","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$year=="2015") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele + don$pop)
res=anova(model)
# Add result to the 'result' file
result_lm_2015_DSDL[num,1:4]=c(i, res$`Pr(>F)`[1:2],summary(model)$r.squared)
}
}
# wrong class
result_lm_2015_DSDL[,2:4]=apply(result_lm_2015_DSDL[,2:4] , 2 , as.numeric)
dim(result_lm_2015_DSDL)## [1] 7341 4
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODs with the genetic map
result_lm_2015_DSDL=merge(map,result_lm_2015_DSDL, by.x=1 , by.y=1, all.y=T)
result_lm_2015_DSDL=result_lm_2015_DSDL[order(result_lm_2015_DSDL$chr, result_lm_2015_DSDL$dist) , ]
# And plot it
plot(-log10(result_lm_2015_DSDL$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2015_DSDL$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2015_DSDL))
num=aggregate(num, by=list(result_lm_2015_DSDL$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test7: DS - 2012 AND 2015
QTL detection without QTL rel for DS only in 2012 only
# Initialize the result table
result_lm_2012_2015_DS=data.frame(matrix(0,0,4))
colnames(result_lm_2012_2015_DS)=c("marker","pval-marker","pval-year","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$pop=="TT") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele + don$year)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_2015_DS[num,1:4]=c(i, res$`Pr(>F)`[1:2],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_2015_DS[,2:4]=apply(result_lm_2012_2015_DS[,2:4] , 2 , as.numeric)
dim(result_lm_2012_2015_DS)## [1] 3544 4
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODs with the genetic map
result_lm_2012_2015_DS=merge(map,result_lm_2012_2015_DS, by.x=1 , by.y=1, all.y=T)
result_lm_2012_2015_DS=result_lm_2012_2015_DS[order(result_lm_2012_2015_DS$chr, result_lm_2012_2015_DS$dist) , ]
# And plot it
plot(-log10(result_lm_2012_2015_DS$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_2015_DS$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_2015_DS))
num=aggregate(num, by=list(result_lm_2012_2015_DS$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test8: DL - 2012 AND 2015
QTL detection without QTL rel for DL only in 2012 only
# Initialize the result table
result_lm_2012_2015_DL=data.frame(matrix(0,0,4))
colnames(result_lm_2012_2015_DL)=c("marker","pval-marker","pval-year","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don=don[which(don$pop=="BX") , ]
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele + don$year)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_2015_DL[num,1:4]=c(i, res$`Pr(>F)`[1:2],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_2015_DL[,2:4]=apply(result_lm_2012_2015_DL[,2:4] , 2 , as.numeric)
result_lm_2012_2015_DL=na.omit(result_lm_2012_2015_DL)
dim(result_lm_2012_2015_DL)## [1] 5851 4
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODL with the genetic map
result_lm_2012_2015_DL=merge(map,result_lm_2012_2015_DL, by.x=1 , by.y=1, all.y=T)
result_lm_2012_2015_DL=result_lm_2012_2015_DL[order(result_lm_2012_2015_DL$chr, result_lm_2012_2015_DL$dist) , ]
# And plot it
plot(-log10(result_lm_2012_2015_DL$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_2015_DL$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_2015_DL))
num=aggregate(num, by=list(result_lm_2012_2015_DL$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test9: DS AND DL - 2012 AND 2015
QTL detection without QTL rel for DS and DL in 2012 and 2015
# Initialize the result table
result_lm_2012_2015_DSDL=data.frame(matrix(0,0,5))
colnames(result_lm_2012_2015_DSDL)=c("marker","pval-marker","pval-year","pval-pop","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele + don$year +don$pop)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_2015_DSDL[num,1:5]=c(i, res$`Pr(>F)`[1:3],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_2015_DSDL[,2:5]=apply(result_lm_2012_2015_DSDL[,2:5] , 2 , as.numeric)
result_lm_2012_2015_DSDL=na.omit(result_lm_2012_2015_DSDL)
dim(result_lm_2012_2015_DSDL)## [1] 7341 5
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODL with the genetic map
result_lm_2012_2015_DSDL=merge(map,result_lm_2012_2015_DSDL, by.x=1 , by.y=1, all.y=T)
result_lm_2012_2015_DSDL=result_lm_2012_2015_DSDL[order(result_lm_2012_2015_DSDL$chr, result_lm_2012_2015_DSDL$dist) , ]
# And plot it
plot(-log10(result_lm_2012_2015_DSDL$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_2015_DSDL$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_2015_DSDL))
num=aggregate(num, by=list(result_lm_2012_2015_DSDL$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
test10: DS AND DL - 2012 AND 2015 - INTERACTION
QTL detection without QTL rel for DS and DL in 2012 and 2015. And I add interaction between marker and year
# Initialize the result table
result_lm_2012_2015_DSDL_inter=data.frame(matrix(0,0,6))
colnames(result_lm_2012_2015_DSDL_inter)=c("marker","pval-marker","pval-year","pval-pop","pval-inter","R2")
num=0
# Run a loop on every markers
for(i in colnames(genotype)){
# build data frame
num=num+1
don=cbind(Y, genotype[,which(colnames(genotype)%in%i)])
colnames(don)[4]="allele"
don$allele=droplevels(don$allele)
# Run the linear model
if(nlevels(don$allele)==2){
model=lm(don$Elisa_blup_AR1 ~ don$allele + don$year +don$pop + don$allele*don$year)
res=anova(model)
# Add result to the 'result' file
result_lm_2012_2015_DSDL_inter[num,1:6]=c(i, res$`Pr(>F)`[1:4],summary(model)$r.squared)
}
}
# wrong class
result_lm_2012_2015_DSDL_inter[,2:6]=apply(result_lm_2012_2015_DSDL_inter[,2:6] , 2 , as.numeric)
result_lm_2012_2015_DSDL_inter=na.omit(result_lm_2012_2015_DSDL_inter)
dim(result_lm_2012_2015_DSDL_inter)## [1] 7341 6
We are supposed to find more or less the same result than with QTL-Rel. Is it true? The only difference is that we use lm and not QTL-Rel, thus we do not take into account the kinship matrix. Let’s check the manathan plot.
# Merge LODL with the genetic map
result_lm_2012_2015_DSDL_inter=merge(map,result_lm_2012_2015_DSDL_inter, by.x=1 , by.y=1, all.y=T)
result_lm_2012_2015_DSDL_inter=result_lm_2012_2015_DSDL_inter[order(result_lm_2012_2015_DSDL_inter$chr, result_lm_2012_2015_DSDL_inter$dist) , ]
# And plot it
plot(-log10(result_lm_2012_2015_DSDL_inter$`pval-marker`) , pch=20 , col=as.numeric(result_lm_2012_2015_DSDL_inter$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" )
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(result_lm_2012_2015_DSDL_inter))
num=aggregate(num, by=list(result_lm_2012_2015_DSDL_inter$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1]) //: =========================================================================================================================
RECAP
# Function to plot
manat_plot=function(data, my_xlab){
plot(-log10(data$`pval-marker`) , pch=20 , col=as.numeric(data$chr) , cex=1.3, xaxt="n", ylab="LOD - scores" , xlab=my_xlab)
abline(h=3.6, col="grey", lwd=1.5)
num=seq(1,nrow(data))
num=aggregate(num, by=list(data$chr), mean , na.rm=T)
axis(1, at=num[,2], labels=num[,1])
}Detection by year and pop separately
par(mfrow=c(2,2))
manat_plot(result_lm_2012_DS, "DS 2012")
manat_plot(result_lm_2012_DL, "DL 2012")
manat_plot(result_lm_2015_DS, "DS 2015")
manat_plot(result_lm_2015_DL, "DL 2015")Grouping population
par(mfrow=c(2,1))
manat_plot(result_lm_2012_DSDL, "DS & DL in 2012")
manat_plot(result_lm_2015_DSDL, "DS & DL in 2015") ## Grouping Years
par(mfrow=c(1,2))
manat_plot(result_lm_2012_2015_DS, "DS in 2012 + 2015")
manat_plot(result_lm_2012_2015_DL, "DL in 2012 + 2015")Grouping Years AND pops
Without interaction
par(mfrow=c(1,1))
manat_plot(result_lm_2012_2015_DSDL, "DS +DL in 2012 + 2015")With interaction
par(mfrow=c(1,1))
manat_plot(result_lm_2012_2015_DSDL_inter, "DS +DL in 2012 + 2015 AND Interaction")